function [M,PI]=mdrkw(B,CF,Q,RHO,LE,MU,reord,nd,lpd,nx,bcrit,schur,toli);

% This program takes the reduced system produced by redkw.m and 
% the decomposition of its dynamic subsystem undertaken by dynkw.m
% and computes the RE solution. Computation can take place using either
% the eigenvalue-eigenvector decomposition or the Schur decomposition
v=version;

if (eval(v(1))>=5)
   npd=length(lpd);
else
   npd=max(size(lpd));   % npd is number of predetermined variables
end

nus=rows(LE); % number of unstable roots.


[ind back]=sort(reord); % reverse transformation

if (nargin<13)
   toli=10^-6;
end

if (nus>0)
   
ny=rows(B);
nf=ny-nd;

% Consider the equation: E d(t+1|t) = W d(t) + PSId(F) E x(t|t)

% The solution strategy is then as follows:  We know that 
%  (1)  LE*Ed(t+1|t) = LE*W d(t) + LE*PSId(F) Ex(t|t)
% We know further that LE*W = MU * LE from the definition
% of left eigenvectors (this holds for any set of eigenvectors).
% Hence, equation (1) is the unstable canonical variables 
% equation.  We can hence solve it "forward" as
%  (2)  LE*d(t) = inv(MU)*LE Ed(t+1|t) -inv(MU)*LE*PSId(F) Ex(t|t)

nlead=cols(CF)/nx;
nlead=nlead-1;

% Calculate CF0*Q*I+CF1*Q*RHO+CF2*Q*RHO^2 + ...CFn*Q*RHO^nlead
% i.e., calculate C(F) Ex(t|t)
R=eye(cols(RHO));
PHI=CF(:,1:nx)*Q*R;
j=0;
while (j<nlead)
 j=j+1;
 R=R*RHO;
 PHI=PHI+CF(:,j*nx+1:(j+1)*nx)*Q*R;
end


% We now have a system in the form 
%   A* E y(t+1)|I(t) = B y(t) + PHI drv(t),
% with the A and B matrices having been "reduced" by dynkw.m
% but we still need to rule out explosive paths, etc. Part of this system
% is E d(t+1)|I(t) = W  d(t) +   PHId drv(t).  Multiplying this by LE we get
%    Eu(t+1)|I(t) =  MU u(t) + LE*PHId drv(t), which are "separated"
% dynamic equations for the nonpredetermined/unstable canonical
% variables u(t) since MU is upper triangular. Set q = LE*PHI2.

% To impose stability, we seek a matrix theta such that
% theta*RHO = MU * theta + q.  There are then two cases
% If MU is diagonal, then the ith row of theta satisfies:
% thetai * RHO = MU(i,i) * thetai + qi => thetai = qi/(RHO-MU(i,i)*eye(RHO);
% If MU is lower triangular, then we must add in some 
% additional terms that involve MU(i,j)*thetaj for j
% =1,2,...i-1.

q = LE*PHI(nf+1:ny,:);

theta=zeros(nus,cols(PHI));
IRHO=eye(cols(RHO));

for i=1:nus;
 theta(i,:)=q(i,:);
   if schur==1
     for j=1:i-1;
       theta(i,:)=theta(i,:)+MU(i,j)*theta(j,:);
     end
   end
  theta(i,:)=theta(i,:)/(RHO-MU(i,i)*IRHO);
end

% We now know that u(t) = LE * d(t) = theta * drv(t).  We solve for lam(t),
% the first nus elements of d(t), as lam(t) = -(Ll\Lk) k(t) + (Ll\theta) drv(t):

if (npd>0)
 lamk=-LE(:,1:nus)\LE(:,nus+1:nd);  
end

lamdr= LE(:,1:nus)\theta;

BT=B;

% We want to impose this solution on the dynamic system, i.e.,
% "substitute out" for lam(t) using the above rule, which we write
% as lam(t) - lamk * k(t) - lamdr * drv(t).  We have equations of the form:
%  B2 * lam(t) +         B3 * k(t) +   PHI           * drv(t) = 0.  Adding  
% -B2 * lam(t) +    B2*lamk * k(t) +   B2*lamdr      * drv(t) = 0,  we get
%  0  * lam(t) + B3+B2*lamk * k(t) +  (PHI+B2*lamdr) * drv(t) = 0

B2 = BT(:,nf+1:nf+nus);
B3 = BT(:,nf+nus+1:ny);
% incorporate the influence of k(t) via lam(t) into system.
if (npd>0)
 BT(:,nf+nus+1:ny)=B3+B2*lamk;
end
% incorporate the influence of drv(t) via lam(t) into system.
PHI=PHI+B2*lamdr;

% patch up the lam(t) influences and the lam(t) equations.
BT(:,nf+1:nf+nus)=zeros(rows(B2),cols(B2));
BT(nf+1:nf+nus,nf+1:nf+nus)=eye(nus);
if (npd>0)
 BT(nf+1:nf+nus,nf+nus+1:nf+nd)=-lamk;
end
PHI(nf+1:nf+nus,:) = -lamdr;

% solve for state space system:
% z(t) =   PI*s(t)
% s(t+1) = M s(t) + e(t+1)
% with s(t) = |k(t)  |   and z(t) = |y(t)|
%             |drv(t)|              |k(t)|
%                                   |x(t)|

ndrv=cols(Q);

PI= [-BT(1:nf+nus,nf+nus+1:ny) -PHI(1:nf+nus,:)
      eye(npd)                  zeros(npd,ndrv)];

% This is a PI matrix for the new ordering, we now create one
% for the old ordering and add exogenous variables. The index
% "back" reverses the earlier ordering.

PI=  [PI(back,:)
      zeros(nx,npd)              Q           ];

M = [BT(nf+nus+1:ny,nf+nus+1:ny)  PHI(nf+nus+1:ny,:)
     zeros(ndrv,npd)              RHO];


else % there is a system with no unstable roots

ny=rows(B);
nf=ny-nd;

% Consider the equation: E d(t+1|t) = W d(t) + PSId(F) Ex(t|t)

% The solution strategy is then very simple because there
% are no variables to be solved "forward". W is simply the 
% feedback matrix and the PSI(F) polynomial governs
% the response to driving variables.


nlead=cols(CF)/nx;
nlead=nlead-1;

% Calculate CF0*Q*I+CF1*Q*RHO+CF2*Q*RHO^2 + ...CFn*Q*RHO^nlead
% i.e., calculate C(F) Ex(t|t)

R=eye(size(RHO));
PHI=CF(:,1:nx)*Q*R;
j=0;

while (j<nlead)
 j=j+1;
 R=R*RHO;
 PHI=PHI+CF(:,j*nx+1:(j+1)*nx)*Q*R;
end

% We now have a system in the form 
%   E y(t+1|t) = B y(t) + PHI drv(t)
% The first ny-np of these equations are for flows
% and are in the form I f(t) + V k(t) + PHIf drv(t)
% and the last np are in the form: k(t+1)=W k(t) +PHId drv(t)


% solve for state space system:
% x(t) =   PI*s(t)
% s(t+1) = M s(t) + e(t+1)
% with s(t) = |k(t)  |   and z(t) = |y(t)|
%             |drv(t)|              |k(t)|
%                                   |x(t)|

ndrv=cols(Q);

PI= [-B(1:nf,nf+1:ny) -PHI(1:nf,:)
      eye(npd)         zeros(npd,ndrv)];

% This is a PI matrix for the new ordering, we now create one
% for the old ordering and add exogenous variables. The index
% "back" reverses the earlier ordering.

PI=  [PI(back,:)
      zeros(nx,npd)              Q           ];

M = [B(nf+1:ny,nf+1:ny)   PHI(nf+1:ny,:)
     zeros(ndrv,npd)                 RHO];


end

% test for the presence of larger imaginary components

IPI=imag(PI);
t=max(max(abs(IPI)));

if t>toli;
   disp(['imaginary component of PI is larger than ', num2str(toli)])
   disp('This tolerance can be adjusted by setting the last argument')
   disp('of the mdrkw.m function, as called in resolkw.m')
   disp('   strike any key to continue')
   pause
end

IM=imag(M);
t=max(max(abs(IM)));

if t>toli;
   disp(['imaginary component of M is larger than ', num2str(toli)])
   disp('This tolerance can be adjusted by setting the last argument')
   disp('of the mdrkw.m function, as called in resolkw.m')
   disp('   strike any key to continue')
   pause
end
M=real(M);
PI=real(PI);

